#################################
# ANALYSIS: TPIs and current constitutions
# Correlations of colonizing power with the constitutionalization of TPIs
# 
# Part of replication of:
# Traditional Institutions in Africa, Past and Present
#
# Political Science Research and Methods
#
# By Clara Neupert-Wentz and Carl Müller-Crepon, 2023
#
##########################################

# PREPARE ANALYSIS ##############

## Define main data
data <- readRDS("data/main_data.rds")

# Load constitutions data from Holzinger et al.
const.df <- readRDS(file.path("data/constitution_df.rds"))


# EFFECT OF COLONIZER ON CONSTITUTIONALIZATION

# Prepare data
modal <- function(x){names(table(x))[which.max(table(x))[1]]}
est.df <- join(const.df,
                 aggregate.data.frame(data[,c( "colonizer.all"), drop = F],
                                      data[,c("iso3c"), drop = F],
                                      FUN = modal),
                 type = "left", by = "iso3c")

# All outcomes
all.out <- c("acknowledgement", "regulation",  "positiveintegration")
all.out.labs <- c("Acknowledgement", "Regulation",  "Integration")

# Center and rescale variables to 0-1
est.df[, all.out] <- t((-1*(apply(est.df[, all.out], 2, min, na.rm = T)) + t(est.df[, all.out])))
est.df[, all.out] <-  t(t(est.df[, all.out]) / (apply(est.df[, all.out], 2, max, na.rm = T)))

# Estimate
model.ls <- lapply(all.out, function(dv){
  felm(as.formula(paste(dv, "~ colonizer.all | 0 | 0 | 0")),
       data = est.df)
})

# Plot coefs
coef.ls <- extract_coef(model.ls)

# Prepare plot
plot.df <- do.call(rbind, lapply(model.ls, function(c){
  data.frame(outcome = rep(c$lhs, nrow(c$beta)),
             treat = gsub("colonizer.all","",rownames(c$beta), fixed = T),
             coef = c$beta[, 1],
             se = diag(c$robustvcv)^.5, stringsAsFactors = F)
}))
plot.df <- plot.df[plot.df$treat != "(Intercept)",]
plot.df$outcome_label <- factor(rep(all.out.labs, each = 4), 
                                levels = all.out.labs, ordered = T)
plot.df$pos <- as.numeric(plot.df$outcome_label )

# Plot
g <- ggplot(plot.df, aes(y = coef, x = pos)) +
  geom_point() +
  geom_errorbar(aes(ymin = coef - 1.96*se,
                    ymax = coef + 1.96*se),
                width=0) +
  geom_hline(yintercept = 0, color = "darkgrey", lty = 2) +
  scale_x_continuous(breaks=c(1:length(all.out)), labels = all.out.labs, limits = c(0.5,3.5))  +
  xlab(NULL) + ylab("Effect of Colonizer on Constitutional Variables") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.35, hjust = 1),
        legend.position = "none") +
  facet_wrap(~ treat, nrow = 1)


png(file.path(fig.path, paste0( "const_colonizer.png")), width = 8, height = 4, unit = "in", res = 400)
par(mar = c(0,0,0,0))
print(g)
dev.off()

